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CONVERSION  FACTORS,  U.  S.  CUSTOMARY  TO  METRIC  (SI) 

UNITS  OF  MEASUREMENT 

U.  S.  customary  units  of  measurement  used  in  this  report  can  he  con¬ 
verted  to  metric  (SI)  units  as  follow: 

_ Multiply _  By  _ To  Obtain _ 

cubic  feet  per  second 
feet 

feet  per  second 
feet  per  second  per  second 
miles  (U.  S.  statute) 
square  feet 


0.02831685  cubic  metres  per  second 

0 . 30U8  metres 

0.30U8  metres  per  second. 

0.30U8  metres  per  second  per  second 

1.6093UU  kilometres 

0 . 0929030U  square  metres 


SAVANNAH  HARBOR  INVESTIGATION  AND  MODEL  STUDY 
REANALYSIS  OF  FRESHWATER  CONTROL  PLAN 


Numerical  Model  Study 

PART  I :  INTRODUCTION 
Description  and  Background 

1.  Savannah  Harbor  is  located  on  the  South  Atlantic  coast,  partly 
in  South  Carolina  and  partly  in  Georgia.  The  harbor  is  formed  by  the 
lower  31  miles*  of  the  Savannah  River,  including  a  bar  channel  extending 
to  the  36-ft  contour  in  the  Atlantic  Ocean.  The  Savannah  River's  head¬ 
water  tributaries  rise  in  the  Blue  Ridge  Mountains  of  North  Carolina 
and  unite  to  form  the  main  stream  which  flows  297  miles  to  the  city  of 
Savannah.  The  Savannah  River  is  tidal  for  approximately  50  miles,  from 
its  mouth  to  Ebenezer  Landing  (Plate  l). 

2.  A  deep-draft  navigation  project  extends  up  Front  River  to  a 
point  just  downstream  from  the  Coastal  Highway  Bridge.  Maintenance  of 
project  depths  has  required  considerable  dredging,  with  the  greater 
portion  being  required  in  the  lower  4  miles  of  Front  River  since  this 
reach  encompasses  the  nodal  point  for  predominance  of  bottom  currents 
under  conditions  of  normal  freshwater  discharge.  As  a  result  of  these 
significant  maintenance  problems,  an  extensive  series  of  physical  model 
investigations  was  conducted  during  the  latter  half  of  the  1950's  and 
early  1960's  to  develop  plans  for  reducing  maintenance  dredging  costs 
for  Savannah  Harbor.** 

*  A  table  of  factors  for  converting  U.  S.  customary  units  of  measure¬ 
ment  to  metric  (SI)  units  is  presented  on  page  3. 

**  U.  S.  Army  Engineer  Waterways  Experiment  Station,  CE.  "Savannah 
Harbor  Investigation  and  Model  Study:  Volume  III,  Results  of  Model 
Investigations,"  Technical  Report  No.  2-580,  Vicksburg,  Miss. 

&.  1961.  "Section  1,  Model  Verification  and  Results  of  General 

Studies." 

b.  1961.  "Section  2,  Tests  of  Improvement  Plans." 

c_.  1963.  "Section  3,  Results  of  Supplemental  Tests." 

<1.  H.  J.  Rhodes  and  H.  B.  Simmons.  1965.  "Section  4,  Results  of 
Tests  of  Increased  Channel  Dimensions." 


3.  The  major  elements  of  the  sediment  basin-tide  gate  plan  devel¬ 
oped  in  the  physical  model  investigations  involve  dredging  a  sediment 
basin  at  the  lower  end  of  Back  River,  construction  of  a  tide  gate  in 
Back  River,  dredging  a  connecting  canal  between  Back  River  and  Middle 
River,  and  channel  improvements  through  McCoombs  Cut  and  into  the  up¬ 
stream  portion  of  Back  River  (Plate  2).  The  concept  of  the  plan  is  for 
the  tide  gate  structure  to  be  open  during  the  flood  portion  of  the 
tidal  cycle  and  then  closed  at  high-water  slack.  This  will  cause  major 
flow  distribution  changes  in  the  lower  channel  complex  since  the  tidal 
prism  in  the  Back  and  Middle  River  areas  must  be  drained  through  the 
Front  River  portion  of  the  system.  The  higher  ebb  velocities  resulting 
in  Front  River  will  tend  to  move  the  sediment  on  through  the  critical 
reach  and  concentrate  much  of  the  shoaling  in  the  sediment  basin  at  the 
lower  end  of  Back  River.  This  will  permit  more  efficient  dredging  and 
disposal  of  the  material. 

4.  The  improved  diversion  channels  through  McCoombs  Cut  and  the 
upper  end  of  Little  Back  River  were  proposed  to  ensure  fresh  water  for 
the  Savannah  Wildlife  Refuge  by  maintaining  salinities  at  the  inlet 
structure  in  Little  Back  River  at  near  preproject  levels.  The  proposed 
flow  diversion  channels  and  improvements  are  shown  in  Plate  3.  The  in¬ 
take  to  the  freshwater  supply  canal  system  that  will  feed  the  wildlife 
refuge  and  some  private  hunting  lands  is  sized  to  supply  300  cfs  for 
inundation  of  the  area  during  the  fall.  The  channel  improvements  were 
designed  for  a  freshwater  inflow  of  4000  cfs  from  Savannah  River  into 
the  back  channels  (2500  cfs  into  Little  Back  River  and  1500  cfs  into 
Middle  River)  assuming  steady-flow  conditions  and  known  channel  bottom 
slopes.  The  impact  on  salinities  in  Little  Back  River  was  checked  in 
the  physical  model  study  by  pumping  the  design  discharges  from  the 
Savannah  River  into  Little  Back  and  Middle  Rivers  (i.e.,  actual  dynamic 
tests  of  the  proposed  channel  enlargements  were  not  conducted  in  the 
physical  model).  For  the  4000-cfs  freshwater  diversion,  salinities 
near  the  freshwater  intake  for  the  wildlife  refuge  were  not  signifi¬ 
cantly  altered  by  the  sediment  trap  plan. 
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The  Problem 


5.  During  1972,  as  the  U.  S.  Army  Engineer  District,  Savannah, 
developed  detail  plans  for  a  new  freshwater  intake  at  the  wildlife 
refuge,  questions  were  raised  concerning  the  effectiveness  of  the  pro¬ 
posed  channel  improvements  in  providing  the  required  fresh  water  in 
Little  Back  River.  A  reevaluation  of  the  complex  unsteady-flow  situation 
that  exists  in  the  upstream  end  of  the  channel  network  near  McCoombs  Cut 
caused  concern  that  net  flows  of  the  design  magnitudes  would  not  be 
diverted  from  Savannah  River  into  the  back  channels.  At  a  conference 
involving  personnel  from  the  Savannah  District;  U.  S.  Army  Engineer 
Division,  South  Atlantic;  Office,  Chief  of  Engineers;  and  the  U.  S.  Army 
Engineer  Waterways  Experiment  Station  (WES),  it  was  decided  that  addi¬ 
tional  studies  of  this  situation  should  be  made  giving  full  consideration 
to  the  unsteady-flow  conditions  which  exist  in  the  upper  channel  complex 
and  their  impact  on  salinities  at  the  location  of  the  proposed  fresh¬ 
water  intake. 


Purpose  and  Scope  of  Study 

6.  The  study  reported  herein  was  undertaken  to  provide  more  reli¬ 
able  information  concerning  flow  and  salinity  changes  in  the  Little  Back 
River  resulting  from  implementation  of  the  sediment  basin-tide  gate 
plan.  In  order  to  provide  this  information,  numerical  models  were  de¬ 
veloped  to  (a)  calculate  unsteady  flows  in  a  complex  network  of  channels 
and  (b)  estimate  time  varying  salinities  throughout  the  system.  Follow¬ 
ing  model  development,  tests  were  run  to  determine  the  change  in  flows 
and  salinities  resulting  from  various  stages  of  development  of  the  sedi¬ 
ment  basin-tide  gate  plan. 

7.  Since  the  Savannah  District  was  almost  ready  to  let  a  contract 
for  construction  of  a  new  freshwater  intake  at  the  wildlife  refuge ,  pre¬ 
liminary  results  from  this  study  were  summarized  in  an  interim  report 
submitted  to  the  district  in  August  1973.  This  report  documents  more 
fully  the  models  developed,  tests  conducted,  and  results  obtained. 
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PART  II:  THE  NUMERICAL  MODELS 


Hydrodynamic  Model 


8.  The  portion  of  the  Savannah  Estuary  of  interest  in  this  study 
includes  a  network  of  channels  as  shown  in  Plates  1  and  2.  Some  ver¬ 
tical  stratification  existed,  in  the  channel  system,  primarily  in  the 
navigation  channel  reach  of  Front  River.  However,  this  study  was  con¬ 
ducted  on  the  assumption  that  a  one-dimensional  calculation  (i.e.  well- 
mixed  system)  would  produce  reliable  indications  of  the  impact  of  the 
plan.  An  unsteady- flow  model  capable  of  treating  the  Savannah  Estuary 
Channel  network  (or  other  channel  networks)  was  developed  for  this  study. 
The  one-dimensional  unsteady-flow  equations  that  formed  the  basis  for 
this  model  are 


and 


3_ 

3x 


uQ  +  gA 


S  2 
AC  R 


0 


(1) 


3A  3& 
3t  3x 


0 


(2) 


and  represent,  respectively,  the  conservation  of  longitudinal  momentum 

and  mass.  In  these  equations 

Q  =  discharge,  cfs 

t  -  time,  sec 

u  =  mean  velocity,  fps 

x  =  distance  along  channel,  ft 

2 

A  =  cross-sectional  area,  ft 

2 

g  =  acceleration  due  to  gravity.,  ft/sec 
n  =  water-surface  elevation,  ft 
C  =  Chezy  coefficient 
R  =  hydraulic  radius,  ft 

Along  the  interior  portions  of  the  various  branches  of  the  channel 
network  these  equations  were  integrated  numerically  using  a 


predictor-corrector  technique.  At  the  boundary  points  of  the  system  and 
at  the  junctions  of  various  branches,  special  treatment  based  on  the 
conservation  of  mass  was  required.  A  more  detailed  discussion  of  the 
hydrodynamic  model  is  presented  in  Appendix  A. 

Salinity  Model 

9.  The  one-dimensional  equation  for  conservation  of  a  conservative 
constituent  is 


8s  +  3s  _  8  „  8s 

8t  U  8x  8x  8x 


(3) 


where 

s  =  conservative  constituent 
u  =  longitudinal  velocity,  fps 
E  =  dispersion  coefficient 

This  equation  was  the  basis  of  the  salinity  model.  Again,  the  equation 
was  numerically  integrated  with  special  treatment  required  at  the  junc¬ 
tions  and  boundary  points.  The  time-varying  flow  field  computed  with 
the  hydrodynamic  model  was  used  as  input  data  to  the  salinity  model. 
More  detailed  discussion  of  the  rationale  for  the  salinity  model  is 
presented  in  Appendix  B. 

Schematization  of  Channel  Network 


10.  The  Savannah  Estuary  Channel  network  between  Fort  Jackson  just 
upstream  from  Elba  Island  and  a  point  near  Ebenezer  Landing  (Plate  l) 
was  simulated  in  the  model.  The  upstream  limit  of  the  model  was  beyond 
the  point  of  tidal  influence.  The  schematic  network  for  preproject 
(1950)  conditions  is  shown  in  Plate  U  and  a  similar  network  for  planned 
conditions  (sediment  basin-tide  gate  plan  fully  implemented)  is  shown 
in  Plate  5.  Table  1  gives  the  lengths  of  each  branch  in  the  networks 
and  the  spatial  increments  (Ax)  used  in  the  numerical  integration  pro¬ 
cess.  In  addition  to  nodes  at  all  channel  junctions,  nodes  are  used 
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at  the  upstream  end  of  the  navigation  channel  because  of  the  large 
change  in  cross  section  and  at  the  tide  gate  location  because  of  the 
necessity  for  gate  opening  and  closing  simulation. 

11.  The  channel  cross-section  data  were  developed  from  hydrographic 
charts  available  from  the  Savannah  District  for  1950  conditions.  Channel 
modifications  since  1950  have  been  substantial,  and  updated  sections 
were  obtained  in  most  of  the  channel  network  to  form  another  set  of 
channel  sections  compatible  with  1972  conditions.  This  data  base  re¬ 
flects  the  sediment  trap  and  the  new  drainage  canal  connecting  Back 
and  Front  Rivers  as  well  as  the  deepening  and  widening  of  the  Front  River 
navigation  channel.  Three  typical  channel  cross  sections  are  presented 
in  Plate  6  to  illustrate  relative  channel  depths  and  widths.  Both  cross- 
section  data  bases  consisted  of  a  series  of  tables  of  channel  widths 
versus  elevations  at  10  or  less  equal  elevation  intervals  plus  the 
elevation  of  the  channel  bottom.  Plate  6  shows  how  a  Front  River  cross 
section  was  divided  into  equidepth  segments. 

Model  Calibration  and  Verification 


12.  Calibration  of  the  hydrodynamic  model  was  based  on  field  data 
collected  in  July  1950.  This  is  some  of  the  data  used  in  verifying  the 
Savannah  Harbor  Physical  Model*  and  was  selected  because  it  approximates 
the  mean  tide  for  the  Savannah  Estuary.  The  calibration  runs  were  made 
using  the  measured  tidal  elevation  time  history  at  Fort  Jackson  as  the 
downstream  boundary  condition  and  the  flow  near  Ebenezer  Landing  as  the 
upstream  boundary  conditions.  Computed  water-surface  time  histories  were 
compared  with  measured  elevations  at  the  tide  gage  locations  shown  in 
Plate  Agreement  was  very  good  with  elevations  differing  by  less  than 
0.5  ft  at  all  stations  and  with  timing  of  high  and  low  water  also  agree¬ 
ing  closely.  Computed  and  measured  elevations  at  Bird  Refuge  on  Back 
River  and  the  ACL  Bridge  in  the  Savannah  River  are  presented  in  Plate  7 
as  examples.  An  effort  was  made  to  verify  some  tide  height  data  obtained 


*  See  footnote  on  page  4. 


in  September  1972  (after  Back  River  was  blocked  with  the  cofferdam  at  the 
the  location  of  the  proposed  gate ) .  Agreement  was  not  as  good  as  ex¬ 
pected,  but  it  was  later  confirmed  by  the  Savannah  District  that  the  gage 
zero  at  Fort  Jackson  (the  downstream  boundary  for  the  numerical  model) 
had  shifted  by  about  0.2  ft.  This  is  felt  to  be  the  primary  reason  for 
poor  verification  of  the  1972  data.  Other  test  runs  were  conducted  of 
the  tidal  hydraulic  model  to  ensure  that  the  model  was  producing  tidal 
elevation  and  velocity  results  consistent  with  prototype  and  physical 
model  results.  It  was  concluded  after  the  several  numerical  experiments 
and  verification  tests  that  the  mathematical  model  satisfactorily  repro¬ 
duced  the  system  and  produced  consistent  and  reliable  results. 

13.  Calibration  of  the  salinity  model  was  more  difficult  because 
of  complexities  related  to  numerically  modeling  the  diffusion  process 
and  the  relatively  limited  data  available  in  the  back  channels.  After 
model  development  and  adjustment,  model  computations  for  preproject 
conditions  with  mean  tide  and  normal  freshwater  inflows  (about  7000  cfs) 
were  in  good  agreement  with  the  measured  high-water  slack  salinity  pro¬ 
file  in  Front  River.  In  the  calibration  runs,  the  dispersion  coefficient 
E  was  set  to  zero  since  the  numerical  dispersion  introduced  by  the 
numerical  integration  procedure  resulted  in  satisfactory  agreement 
between  computed  and  measured  salinities.  Checks  with  other  available 
data  also  appeared  reasonable  and  led  to  confidence  that  model  calcula¬ 
tions  would  provide  reasonable  estimates  of  salinities  throughout  the 
system  for  different  stages  of  plan  implementation  and  with  the  fully 
implemented  plan.  A  more  extensive  salinity  verification  was  not  con¬ 
sidered  practical  due  to  time  limitations  and  the  requirement  for  exten¬ 
sive  field  data  in  the  channel  network  simulated  by  the  model. 
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PART  III:  TEST  PROGRAM  AND  RESULTS 


Test  Conditions 


14.  The  model  test  program  was  designed  to  provide  estimates  of 
flow  and  salinity  conditions  at  different  stages  of  implementation  of 
the  sediment  basin-tide  gate  plan.  Test  1  represented  preproject  condi¬ 
tions  and  Test  2  represented  the  situation  that  existed  at  the  time  of 
the  study  (cofferdam  blocking  Back  River  at  tide  gate  location  and  a 
connecting  canal  dredged  between  Back  River  and  Middle  River).  Test  3 
simulated  the  situation  that  would  exist  with  the  tide  gate  completed 
and  operating  without  any  channel  improvements  to  McCoombs  Cut  and  the 
upstream  reach  of  Little  Back  River.  The  completed  plan  as  proposed 
was  simulated  in  Test  4. 

15.  The  four  test  conditions  cited  above  had  been  run  at  the  time 
of  the  preliminary  report.  Subsequent  to  that  report,  two  additional 
tests  were  run  to  verify  that  net  downstream  flows  in  Little  Back  River 
could  be  increased  by  additional  changes  in  the  back  channels.  Test  UA 
was  a  repeat  of  Test  4  except  that  the  new  connecting  canal  between  Back 
River  and  Middle  River  was  dredged  to  el  -31  rather  than  el  -9,  as  used 
in  Test  4.  Test  5  was  run  to  determine  if  blocking  Middle  River  near 
its  upstream  intersection  with  Back  River  would  increase  flow  into 

Back  River.  Each  test  condition  is  briefly  summarized  in  Table  2. 

1 6.  All  tests  were  quasi-steady-state  tests  (i.e.,  repetitive  tide 
and  salinity  at  the  downstream  boundary  with  a  constant  upstream  inflow) 
using  the  same  upstream  and  downstream  boundary  conditions.  The  input 
tide  at  the  downstream  boundary  of  the  system  (Fort  Jackson)  was  the 
mean  tide  used  in  the  earlier  physical  model  tests  (Plate  8).  A  constant 
freshwater  inflow  of  7300  cfs  was  used  as  the  upstream  boundary  condi¬ 
tion  in  the  Savannah  River  near  Ebenezer  Landing.  Union  Creek,  which 
flows  into  McCoombs  Cut,  was  assumed  for  this  study  to  be  a  storage 
basin  that  fills  and  drains  with  a  reflecting  boundary.  This  is  con¬ 
sidered  a  reasonable  assumption  to  provide  comparative  test  results.  In 
the  salinity  model,  the  input  salinity  at  Fort  Jackson  (Plate  12)  was 


taken  from  field  data  near  Elba  Island  on  20  July  1950  and  represents 
an  average  over  depth.  The  tidal  conditions  during  that  period  approxi¬ 
mated  mean  tide  conditions. 

Discussion  of  Results 

Tests  1-4 

17.  The  results  of  Tests  1-4  presented  in  this  report  and  dis¬ 
cussed  in  subsequent  paragraphs  differ  slightly  from  those  presented  in 
the  preliminary  report  furnished  the  Savannah  District  in  August  1973. 
The  differences  result  from  the  correction  of  two  of  the  data  points 
defining  the  mean  tide  curve  at  Fort  Jackson  that  was  used  as  the 
driving  downstream  boundary  conditions  for  all  tests.  This  correction 
smoothes  out  some  local  irregularities  in  the  computed  results  and 
changes  some  of  the  numerical  values  cited  but  does  not  impact  the  con¬ 
clusions  presented  in  the  preliminary  report. 

18.  Computed  tidal  elevations  at  selected  locations  in  the  channel 
complex  are  presented  in  Plates  8  and  9.  The  impact  of  the  different 
test  conditions  on  tidal  elevations  at  Bull  Street  on  Front  River  is 
insignificant.  Farther  upstream  on  the  main  channel  at  ACL  Bridge  (just 
upstream  from  McCoombs  Cut),  the  different  test  conditions  do  result  in 
measurable  changes  in  tidal  range  and  phasing.  However,  the  difference 
between  Test  1  (preproject)  and  Test  4  (plan  fully  implemented)  is 
quite  small.  At  Bird  Refuge  on  Back  River  the  changes  are  more  pro¬ 
nounced.  Closing  of  Back  River  with  the  cofferdam  in  Test  2  and  with 
the  tide  gate  in  Tests  3  and  4  delays  the  ebb  portion  of  the  cycle  by 
about  three  quarters  of  an  hour  at  this  location.  The  phase  shift  re¬ 
sults  from  the  longer  flow  path  required  when  Back  River  is  blocked 
during  ebb.  With  the  plan  fully  implemented  (Test  4),  the  flood  portion 
of  the  tidal  cycle  at  Bird  Refuge  is  advanced  by  about  15  minutes  as  a 
result  of  faster  tidal  wave  propagation  through  the  improved  channels, 
the  plan  also  reduces  the  tidal  range  about  0.6  ft  at  Bird  Refuge. 

19.  Discharges  calculated  through  the  tidal  cycle  at  selected 
stations  are  plotted  in  Plates  10  and  11.  Discharges  through  McCoombs 
Cut  and  at  Bird  Refuge  on  Back  River  (Plate  10 )  are  important  in 
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analyzing  the  effect  of  the  plan  on  flow  and  salinity  conditions  in  Back 
River  and  are  summarized  "below: 


McCoombs  Cut 

Bird  Refuge 

on  Back  River 

Peak  D.S.  Net  D.S. 

Peak  D.S. 

Net  D.S. 

Test 

Discharge 

Discharge 

Discharge 

Discharge 

Condition 

cfs 

cfs 

cfs 

cfs 

1 

530 

135 

4l6o 

190 

2 

1370 

440 

2920 

3 

1140 

435 

2850 

75 

4 

3230 

1350 

3190 

680 

In  Test  1  (preproject  conditions)  the  net  downstream  flow  at  Bird 
Refuge  is  larger  than  the  net  flow  through  McCoombs  Cut  because  there 
is  a  net  upstream  flow  in  Middle  River.  When  compared  with  Test  1  dis¬ 
charges,  Test  3  (before  channel  improvements  in  McCoombs  Cut  and  upper 
Back  River)  results  in  some  increase  in  flow  through  McCoombs  Cut  but 
decreases  the  net  downstream  flow  at  Bird  Refuge  on  Back  River.  Im¬ 
provement  of  the  channels  in  the  McCoombs  Cut  area  (Test  4)  approxi¬ 
mately  triples  the  net  downstream  discharge  through  McCoombs  Cut  and 
increases  the  net  downstream  discharge  at  Bird  Refuge  by  a  factor  of 
about  nine  in  comparison  with  computed  discharges  in  Test  3.  Computed 
discharges  at  Bull  Street  on  Front  River  and  at  cross  tides  on  Back 
River  are  plotted  in  Plate  11  to  document  the  flow  changes  in  other 
parts  of  the  system. 

20.  Plates  12  and  13  present  computed  salinity  variations  through 
the  tidal  cycle  at  selected  stations.  High-water  slack  salinity  pro¬ 
files  in  Back  River  and  Front  River  are  shown  in  Plates  l4  and  15 »  re¬ 
spectively.  Operation  of  the  tide  gate  structure  in  Test  3  increases 
the  peak  salinity  at  Bird  Refuge  on  Back  River  to  about  4.8  ppt.  The 
proposed  improvements  in  the  channel  through  McCoombs  Cut  and  in  the 
upper  end  of  Back  River  (Test  4)  reduce  the  peak  salinity  to  about 
2.2  ppt.  The  computed  peak  salinity  for  preproject  conditions  (Test  l) 
was  0.9  at  this  location.  The  higher  salinities  in  Back  River  with 
tide  gate  operation  result  from  the  trapping  of  more  saline  water  in 
Back  River  when  the  tide  gate  is  closed  at  high-water  slack.  The  net 
effect  of  this  is  to  transfer  the  salinity  at  the  tide  gate  location 


during  the  flood  portion  of  the  cycle  (about  6  ppt)  to  the  location  of 
the  new  dredged  canal  connecting  Back  River  with  Middle  River.  How 
rapidly  the  salinities  decrease  upstream  from  that  point  on  Back  River 
depends  upon  the  magnitude  of  the  freshwater  inflows  into  Back  River. 
Study  of  the  salinity  profile  in  Plate  l4  shows  the  effect  of  the  changes 
in  discharges  in  the  back  channels  discussed  in  the  preceding  paragraph. 

21.  Salinity  changes  in  Front  River  are  less  important  but  are 
worth  noting  (Plate  15).  Implementation  of  the  tide  gate  plan  lowers 
the  high-water  slack  salinity  profile  between  Fort  Jackson  and  Bull 
Street  and  significantly  raises  the  profile  from  Bull  Street  to  a  point 
just  upstream  from  the  end  of  the  navigation  channel.  Peak  salinities 
at  the  Sugar  Refinery  location  just  upstream  from  the  canal  connection 
with  Back  River  are  increased  from  essentially  zero  for  preproject  con¬ 
ditions  to  about  2.5  ppt  for  the  fully  implemented  plan  (Plate  13). 

22.  Based  on  the  results  of  these  four  tests,  as  reported  in  the 
preliminary  report  on  the  study,  it  was  suggested  that  the  Savannah 
District  proceed  with  implementation  of  the  plan  as  originally  designed. 
This  suggestion  was  based  on  the  following  rationale.  While  the  net 
freshwater  inflow  into  Back  River  was  not  as  high  as  that  originally 
estimated  (computed  at  about  680  cfs  compared  with  2500  cfs  used  in  an 
earlier  physical  model  study) ,  the  calculations  indicated  that  the  peak 
salinity  at  Bird  Refuge  would  only  increase  from  about  0.9  ppt  for  pre¬ 
project  conditions  to  approximately  2.2  ppt  with  the  fully  implemented 
plan.  No  specific  guidance  had  been  provided  concerning  the  acceptable 
peak  salinity  at  Bird  Refuge,  but  it  was  felt  that  the  projected  in¬ 
crease  would  probably  not  be  objectionable.  Also,  should  monitoring  of 
the  situation  at  Bird  Refuge  indicate  the  salinity  increase  to  be  detri¬ 
mental,  the  situation  can  be  reversed  to  essentially  preproject  condi¬ 
tions  by  leaving  the  tide  gate  open  throughout  the  tidal  cycle.  This 
would  provide  a  short-term  solution  and  it  was  felt  that  further  changes 
and/or  improvements  in  the  back  channels  could  be  made  to  increase  the 
freshwater  inflow  into  Back  River.  A  short  time  after  the  preliminary 
report  was  prepared,  the  Savannah  District  authorized  WES  to  make  two 
additional  runs  to  verify  that  larger  freshwater  flows  in  Back  River 


could  "be  obtained  if  necessary.  The  results  of  these  tests  are  briefly 
described  in  subsequent  paragraphs. 

Tests  4a  and  5 


23.  A  number  of  possible  channel  modifications  which  offer  poten¬ 
tial  for  increasing  downstream  net  flow  in  Back  River  were  considered: 

a.  Closing  of  Middle  River  in  its  upper  reaches. 

b.  Closing  of  Middle  River  in  its  lower  reaches. 

£.  Construction  of  a  flow  restriction  in  the  upper  reaches  of 
Middle  River. 

d.  Further  extension  of  the  channel  modifications  in  Little 
Back  River  to  Bird  Refuge. 

e_.  Channel  improvements  in  the  connection  between  Front  and 
Back  Rivers. 

£.  Modification  in  gate  operation  to  reduce  salinities  such 
as  gate  closure  for  one  or  two  tidal  cycles  per  week. 

Modifications  e_  and  a  were  selected  for  testing  as  Tests  4a  and  5> 
respectively.  Discharges  calculated  at  McCoombs  Cut  and  Bird  Refuge 
throughout  the  tidal  cycle  are  plotted  in  Plate  l6  along  with  results 
from  Test  4.  Peak  and  net  discharges  are  summarized  below. 


McCoombs  Cut _  Bird  Refuge  on  Back  River 


Test 

Condition 

Peak  D.S. 
Discharge 
cfs 

Net  D.S. 
Discharge 
cfs 

Peak  D.S. 
Discharge 
cfs 

Net  D.S. 
Discharge 
cfs 

4 

3230 

1350 

3190 

680 

4a 

3160 

1540 

3850 

930 

5 

3095 

890 

3460 

835 

Deepening  the  connecting  canal  between  Back  River  and  Front  River  in- 

creased  the 

net  downstream  discharge  in  McCoombs  Cut  and  at  Bird  Refuge 

by  about  15 

and  35  percent. 

respectively. 

Blocking  Middle 

River  near 

the  upstream  end  actually  reduced  net  flow  through  McCoombs  Cut  but 
increased  the  net  downstream  flow  in  Back  River  by  about  20  percent. 

24.  The  impact  of  the  increased  net  flow  on  peak  salinity  at  Bird 
Refuge  is  indicated  in  the  following  tabulation. 

Test  Condition  Peak  Salinity  at  Bird  Refuge,  ppt 

4 

4a 

5 


15 


2.2 

2.3 

2.0 


Blocking  Middle  River  (Test  5)  resulted  in  a  small  decrease  in  peak 
salinity  to  about  2.0  ppt.  However,  deepening  the  connecting  canal 
between  Back  River  and  Middle  River  (Test  kA)  resulted  in  a  slight 
increase  in  peak  salinity  even  though  this  change  produced  larger  in¬ 
creases  in  net  downstream  flow  at  McCoombs  Cut  and  Bird  Refuge  than  the 
change  made  in  Test  5.  The  slight  increase  in  peak  salinity  is  attrib¬ 
uted  to  the  longer  period  of  flood  currents  into  Back  River  which  re¬ 
sulted  with  this  test  condition.  These  channel  modifications  had 
negligible  effects  on  salinities  in  Front  River. 

25.  Results  of  Tests  UA  and  5  verify  that,  if  necessary,  additional 
channel  modifications  can  be  developed  to  increase  freshwater  flows 
into  Back  River.  However,  they  also  show  that  channel  modifications 
must  be  carefully  evaluated  to  ensure  that  the  larger  net  flows  in  Back 
River  also  reduce  peak  salinities  at  Bird  Refuge.  No  attempt  was  made 
to  identify  the  optimum  channel  modification  since  the  need  for  further 
modification  has  not  been  indicated. 
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PART  IV:  CONCLUSIONS  AND  RECOMMENDATIONS 


26.  The  numerical  calculations  reported  herein  suggest  the  follow¬ 
ing  conclusions  and  recommendations. 

a.  The  net  freshwater  inflow  into  Back  River  with  the  sedi¬ 
ment  basin-tide  gate  plan  fully  implemented  will  not  be 
as  large  as  had  earlier  been  estimated.  The  dynamic 
calculations  indicate  a  net  downstream  flow  of  about  680 
cfs  compared  with  2500  cfs  used  in  the  earlier  physical 
model  study. 

b.  Pull  implementation  of  the  plan  will  cause  a  relatively 
small  increase  in  peak  salinities  at  Bird  Refuge  on  Back 
River,  from  about  0.9  ppt  for  preproject  conditions  to 
approximately  2.2  ppt  for  the  fully  implemented  plan. 

c_.  Since  no  specific  guidance  is  available  concerning  accept¬ 
able  peak  salinities,  monitoring  of  salinity  conditions  at 
Bird  Refuge  is  recommended  as  the  plan  is  implemented. 
Should  the  relatively  small  salinity  increase  prove  detri¬ 
mental,  the  situation  can  be  reversed  to  essentially  pre¬ 
project  conditions  by  leaving  the  tide  gate  open  through¬ 
out  the  tidal  cycle.  This  type  of  operation  would  only  be 
appropriate  for  short  periods  of  time  or  until  other 
modifications  could  be  made  since  continuous  operation 
with  the  tide  gate  open  would  make  the  sediment  trap  plan 
ineffective. 

cl.  If  necessary,  the  net  downstream  discharge  in  Back  River 
could  be  increased  by  further  modifications  in  the  back 
channels.  Two  of  the  several  possible  channel  changes 
were  tested  to  verify  that  the  net  downstream  flow  could 
be  increased,  but  no  attempt  was  made  to  identify  the 
optimum  channel  modification  since  the  need  for  further 
modification  has  not  been  indicated. 
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Table  1 

Total  Branch  and  Incremental  Lengths 


Preproject  (1950) 


Conditions 


Branch 

No. 

Length 

ft 

Ax 

ft 

1 

6,700 

1,675 

2 

45,600 

2,682 

3 

86,800 

2,553 

4 

15,100 

1,510 

5 

12,600 

1,400 

6 

8,000 

1,600 

7 

19,800 

2,475 

8 

17,400 

1,933 

9 

1,100 

550 

10 

6,100 

1,220 

11 

14,200 

2,840 

12 

114,655 

U,985 

Present  (l972) 


Conditions 


Branch 

Length 

Ax 

No. 

ft 

ft 

Changed  Branches 


3 

23,200 

2,700 

4 

5,280 

1,760 

5 

7,500 

1,500 

New 

Branches 

13 

12,765 

2,553 

15 

2,600 

1,300 

16 

50,540 

2,660 

17 

5,190 

1,730 

18 

9,810 

1,635 

Proposed  Channel 
Improvements 
Branch  Length  Ax 

No.  ft  ft 


Changed  Branches 

10  It,  800  960 

16  44,460  2,340 


NOTE:  ALL  OTHER  BRANCHES  REMAINED  AS  THE 
1950  CONFIGURATION. 


Table  2 

Test  Conditions 


Test 

No. 

Channel 

Configuration 

Gate 

Diversion 

Channels 

Notes 

1 

1950 

None 

Unimproved 

Preproject 

2 

1972 

Closed 

Unimproved 

Existing 

3 

1972 

Operating 

Unimproved 

No  dredging  . 

H 

1972 

Operating 

Improved 

Planned  conditions 

kA 

1972 

Operating 

Improved 

Connecting  canal 
between  Back 
River  and  Middle 
River  deepened 

5 

1972 

Operating 

Improved 

Middle  River 
blocked  at 
upstream  end 

IMPROVED  CHANNEL 


SAVANNAH  RIVER 


PLATE  3 


SAVANNAH  RIVER 


BRANCH  U2\ 


-UNION  CREEK 


MC  COOMBS 
\CUT ^ 


ACL  BRIDGE 


"’nw 


'MIDDLE  RIVER 


COASTAL  HIGHWAY  BRIDGE 


MIDDLE  RIVER 


'  BIRD  REFUGE 


^-little  back  river 


SUGAR  REFINERY  ( 


KINGS  ISLAND  I 


CROSS  TIDES 


FRONT  RIVER  ^ 


-BACK  RIVER 


BULL  STREET 


BACK  RIVER 


©  TIDE  STATION 
INPUT  TIDE 


FORT  JACKSON 

"^TO  OCEAN 


SAVANNAH  ESTUARY 
NUMERICAL  MODEL 
SCHEMATIC  NETWORK 

PREPROJECT  (1950)  CONDITIONS 


PLATE  4 


SAVANNAH  RIVER 


NUMERICAL  MODEL 
SCHEMATIC  NETWORK 
PLANNED  CONDITIONS 


TYPICAL  CROSS  SECTIONS 


SAVANNAH  ESTUARY 
NUMERICAL  MODEL 

MODEL  ADJUSTMENT 


ELEVATION,  FT- PULASKI  OATUM 
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INPUT  TDE  AT  FT.  JACKSON 


PLATE  8 


TIDAL  COMPUTATIONS 

FT.  JACKSON  AND 
FRONT  RIVER  AT  BULL  STREET 
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PLATE  15 


TIME,  HOURS  AFTER  MOON'S  TRANSIT  OF  OJARANTINE  MERIDIAN 


LEGEND  . 

TEST  NQ  SYMBOL 


TIOE .  MEAN 

FRESHWATER  DISCHARGE  7300  CFS  DISCHARGE  COMPUTATIONS 

MtCOOMBS  CUT  AND 
LITTLE  BACK  RIVER  AT  BIRD  REFUGE 
TESTS  4,  4A,  AND  5 


L 


PLATE  16 


DISCHARGE  1000  CFS 


APPENDIX  A:  THE  HYDRODYNAMIC  MODEL 

Introduction 

1.  The  numerical  study  of  the  upper  reaches  of  the  Savannah 
Estuary  was  based  on  an  application  of  the  one-dimensional  long-wave 
equations.  The  computer  code  which  was  developed  is  also  applicable  to 
the  solution  of  other  problems  where  the  situation  is  topologically 
similar  and  where  long-period  tidal  flows  or  natural  floods  are  of 
interest.  The  program  is  not  suited  to  investigations  in  which  the 
movement  of  bores  must  be  considered. 

2.  Topologically,  the  upper  reaches  of  the  Savannah  Estuary  con¬ 
sist  of  a  network  of  interconnected  channels  and  the  computer  program 
developed  to  support  this  investigation  was  written  to  allow  the  deter¬ 
mination  of  tidal  motion  for  entirely  arbitrary  channel  networks.  As 

a  matter  of  definition,  we  shall  refer  to  a  reach  of  channel  between 
two  junctions  as  a  branch  and  to  the  junctions  as  nodes;  the  extremity 
of  a  channel,  not  at  a  node,  is  referred  to  as  a  boundary.  The  treat¬ 
ment  that  was  applied  to  the  determination  of  flow  conditions  at  a 
junction  is  quite  general  and  the  code,  therefore,  was  written  to  permit 
the  connection  of  two  or  more  branches  at  a  junction  or  a  node.  The 
junction  algorithm  is  also  applicable  to  the  determination  of  tidal 
flows  at  places  where  there  is  a  sudden  change  in  channel  geometry  and 
nodes,  therefore,  were  put  at  locations  where  sudden  changes  in  channel 
geometry  occurred. 

3.  Several  methods  of  deriving  the  one-dimensional  long-wave 
equations  are  available.  The  control  volume  approach  was  particularly 
useful  in  making  decisions  pertinent  to  the  design  of  the  numerical 
model.  The  key  result  of  this  reasoning  being  that  the  long-wave  approx¬ 
imate  momentum  equation  is  not  applicable,  i.e.,  not  descriptive  of  the 
phenomenon,  at  the  junction  of  channels  or  at  locations  where  there  are 
sudden  changes  in  channel  configuration.  This  is  the  case  because  in 
these  instances,  there  are  forces  that  influence  the  flow  but  that  are 
not  accounted  for  in  the  approximation.  This  conclusion 
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is,  of  course,  not  original  (see  Lamb).* 


Numerical  Hydraulics  for  Channel  Networks 
4.  The  one-dimensional  long-wave  equations  used  in  this  study  were 


and 


ifi  =  _ 
at 


la  _  gftlqj 

ac2r 


9A  _  _  IQ 

at  “  ax 


(Al) 


(A2) 


where 

Q  =  discharge 

u  =  mean  velocity,  i.e.,  the  velocity  averaged  over  the  cross 
section 

A  =  cross-sectional  area 
q  =  water-surface  elevation 
C  =  Chezy  coefficient 
R  =  hydraulic  radius 

g  =  gravitational  acceleration  constant 
x  =  longitudinal  distance  along  the  channel 
t  =  time 


In  the  computation  C  was 


n 


(A3) 


*  For  the  case  of  inviscid  tidal  motion  in  a  channel  of  constant  width 
with  a  smooth  bottom,  the  forces  accounted  for  by  the  momentum  equa¬ 
tion  are  the  pressure  forces  acting  at  the  two  ends  of  the  control 
section  and  the  pressure  force  acting  along  the  bottom  of  the  channel. 
If  there  is  a  step  in  the  channel,  an  additional  force  not  included 
in  the  above  acts  on  any  control  reach  lying  across  the  sudden  rise. 
Lamb  treats  this  problem  which  is  that  of  wave  reflection  and  trans¬ 
mission.  (Sir  Horace  Lamb.  1932.  Hydrodynamics ,  6  ed. ,  Dover 
Publications. ) 
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where  n  is  Manning's  roughness  coefficient. 

5.  The  numerical  approximation  of  the  one-dimensional  long-wave 
equations,  which  was  used  in  this  analysis,  may  he  derived  in  several 
ways.  The  one  that  follows  elucidates  the  conservative  characteristics 
of  the  approximations.  Suppose  that  a  branch  of  the  channel  network 
is  divided  into  an  even  number  of  equal  segments,  2n  .  There  are  then 
(including  both  ends)  2n  +  1  stations  along  the  channel.  We  suppose 
these  to  be  located  at  x^  ,  x^  ,  x^  ,  ...  ,  x2n+i  *  *’or  convenience 
we  set 


Qi  =  Q  (xiSt)  ,  Q±  =  f^-  Q(x±,  t) 

•  3 

—  A  (x^,"b)  ,  A(x^,  t ) 


(AU) 


Integrating  Equation  A2  and  applying  Leibnitz's  rule  results  in 


A  dx  =  Qi_1  -  Qi+1  ;  i  =  2,  3, 


2n  -  1  (A5) 


If  the  integral  on  the  left  of  Equation  A5  is  approximated  by  the  mid¬ 
point  rule,  the  result  is 


^Qi-1  "  Qi+1^ 


i=2,  3,  ...  ,  2n  -  1 


where  <5x  is  the  uniform  length  of  the  segments.*  Similarly,  Equa¬ 
tion  A1  becomes 


*  The  use  of  the  midpoint  rule  is  arbitrary  in  the  sense  that  Simpson's 
rule  could  also  be  used,  or  if  a  different  choice  of  elements  were 
made  any  other  numerical  quadrature  formula.  Should  Simpson's  rule 
have  been  used,  the  result  would  have  been 


Z  Vi 


2  • 

+  3  Ai 


"  Qi-1^ 


which  is  two  orders  of  magnitude  more  accurate  than  Equation  AU.  More 
computations  would,  however,  be  required. 


S.  “  26x  [^U<^i-1  "  ^uQ>i+l]  +  26x  ^ni-l  “  ni+l^  ~  (6  .„2. 
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1  (AT) 

i  =  2 ,  3 ,  ...  ,  2n  -  1 


when  the  approximation 


Xi+1  _xi+l 


JfG  dx  =  J F±G  dx  +  0(5x3) 


(A8) 


i-1 


Xi-1 


has  been  applied.  The  analytical,  procedure  just  applied  has  reduced  the 
problem  of  solving  a  pair  of  partial  differential  equations  to  that  of 
solving  a  system  of  ordinary  differential  equations  whose  solution 
approximates  that  of  the  original  problem.  The  procedure  has  numerous 
applications  and  is  referred  to  as  Kantorovich's  method  or  the  method 
of  undetermined  functions. 

6.  The  set  of  double  intervals,  x..  ^  to  xi+1  ,  covers  the 
branch  in  the  sense  that  there  is  no  open  set  along  the  branch  not  lying 
in  an  interval  and  that  every  connected  open  set  along  the  branch  which 
does  not  contain  an  x^  lies  in  a  unique  interval.  Also,  Equation  A6 
clearly  implies  that  the  rate  of  mass  flux  out  of  the  double  interval 

xi  2  ’  xi  e<lua^  the  rate  of  mass  flux  into  the  double  interval 
xi  ’  xi+2  *•  Hence,  it  follows  that  mass  is  conserved,  i.e.,  none  is 
lost  or  gained  by  the  approximation  given  in  Equation  A6. 

3 

7.  Momentum ,  however,  is  conserved  only  up  to  an  order  of  <5x 
because  of  the  approximation  indicated  by  Equation  A8.  It  would  be 
feasible  to  use  a  more  accurate  approximation  than  Equation  A8,  but 
this  would  have  been  pointless  because  the  midpoint  rule  used  in  de- 
riving  Equations  A6  and  AT  is  only  accurate  to  the  order  Sx  . 


Numerical  Integration 


8.  The  2n  -  1  first  order  ordinary  differential  equations  of 


Equation  A6  and  the  2n  -  1  first  order  ordinary  differential  equations 
of  Equation  A7  must  be  solved.  This  could  be  accomplished  with  any  of 
a  number  of  the  various  numerical  methods  for  solving  systems  of  first- 
order  ordinary  differential  equations.  The  two-point  predictor-corrector 
method  was  used  in  this  study.  To  explain  the  method  we  first  note  that 
Equations  A6  and  A7  could  be  written  as  a  system  of  Un  -  2  ordinary 
differential  equations 


z  =  f  (t,z)  (A 9) 

where  z  ,  f  ,  and  z  denote  the  vectors 

z  =  (z1,  z2,  ...  , 

f  =  (fr  f2,  ...  ,  f^)  (A10) 

z  =  (zlf  z2,  ...  ,  z^) 

The  predictor  algorithm  is 

+  {t  (U1) 

where 

Zj  =  z(tj)  ,  fj  =  f(tj,  zj)  (A12) 

and 

«  -  tJW  -  tj  (A13) 


The  corrector  algorithm  is 


where 


r*  =  f(t 

:j+ 1  v 


j+i*  j+i 


To  apply  these  algorithms  z  must  be  known  at  the  time  tj  .  The 
quantity  z*+^  is  the  prediction  of  z  for  the  time  tj+^  .  After  it 
is  computed,  an  improved  or  corrected  estimate  of  z  at  time  tj+1  is 
made  by  applying  Equation  AlU. 

9.  The  predictor-corrector  method  was  employed  because  it  was  known 
that  the  method  is  well  suited  to  the  integration  of  ordinary  differen¬ 
tial  equations  describing  oscillations.  No  comparison  of  the  computa¬ 
tional  efficiency  of  this  algorithm  with  others  available  was  made. 

Our  current  inclination  for  future  studies  would  be  to  make  a  more 
thorough  investigation  of  this  with  the  view  of  using  the  predictor- 
corrector  scheme  for  an  initial  time-step  and  the  central  difference 
scheme 


=  - j— 1  * 


thereinafter. 


The  Junction  of  Channels 


10.  In  preceding  sections,  a  numerical  method  was  given  for  esti¬ 
mating  the  quantities  Q  and  A  at  time  t  +  6t  given  the  values  of 
these  quantities  at  time  t  .  The  numerical  formulas  do  not,  however, 
apply  to  either  the  first  or  last  station  along  a  branch.  Computation 
of  Q  and  A  at  the  ends  of  a  branch  is  considered  in  this  section. 

11.  Q  and  A  may  be  assumed  to  be  known  at  all  stations  at  time 

t  and  using  numerical  algorithms  indicated  in  preceding  sections  Q 
and  A  may  be  computed  at  times  t  +  6t  at  all  interior  stations. 
Hence,  at  a  typical  node  0  as  depicted  in  Figure  Al,  the  discharges 
and  cross-sectional  areas  are  known  as  stations  (a,l)  ,  (b,2)  ,  and 

(c,2)  which  are  interior  stations. 


Figure  A1 

Our  notation  is  that  the  first  variable,  i.e.,  a,  b,or  c,  denotes 
the  branch  and  the  second  variable  denotes  the  station.  It  is  assumed 
that  there  is  no  storage  at  the  node;  and  conservation  of  mass, 
therefore,  requires  that 


Q. 


'a, 2  "  S>,1  " 


Q 


c,l 


=  0 


It  was  also  assumed  that 


(A15) 


V,2  =  \,1  ~  nc,l 


(Al6) 


vhich  is  an  assumption  that  there  was  no  significant  head  loss  at  the 
Junctions. 

12.  These  compatibility  conditions  are  used  along  with  the  integral 
form  of  the  continuity  equation  in  finite  difference  form 

&V-  —  — 

St  *  JQi  -  <"T> 


where 


«  :  5  +  <W>  (Al8) 

that  is,  the  rate  of  change  in  volume  is  equal  to  the  difference  in 
discharges  through  the  node.  The  trapezoidal  rule,  which  is  accurate  to 
third  order,  was  used  for  volume  computations. 
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AT 


or 


fix 

»e  :  (jAlfl  *  «Aa>2) 


(A20) 


For  the  "Y"  shaped  control  volume  of  Figure  Al,  the  continuity  equation 
is  then  written  as 


<Sx 


fix. 


+  SA.,2>  —  +  (8Ab,l  *  SAb,2>  —  4  (SAc,l  *  {Ae,2>  T 
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The  values  of 


fiA 


fix 
_ 8 

a,l  2 
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fix 
_ c 
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are  known  at  time  level  t  +  fit  and  it  is  convenient  to  define  the 
quantity  by 


fix 


fij 


fix 


^=6Aa,l  IT  +6\,2^+6Ac,2- 


(A22) 


Similarly,  the  discharge  values  appearing  in  Equation  A21  are  known 
at  time  level  t  +  fit  and  we  define 


S«  ‘  !Vl  -  V  -  «0,2>  5t 


(A23) 


Thus,  Equation  A21  reduces  to 
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Because  of  Equation  Al6  we  may  define  a  function  F^  for  each  node 
as  a  function  of  water-surface  elevation  n  , 

fix  fix.  fix 

FB^  =  Aa,2  ~~2~  +  ^,1  ~2~  *  Ac,l  ~2~ 


(A25) 
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The  functions  F0  are  known  functions,  i.e.,  data  describing  the  varia- 
p 

tion  of  cross-sectional  area  with  depth  are  data  input  into  the  computa¬ 
tion.  Introduction  of  this  definition  reduces  Equation  A2U  to 


SP6(n)  =  S,  -  S, 


where 


(A26) 


Fg[n(t  +  St )]  =  FB[n(t)]  +  6Fg(n)  (A27) 

Thus,  the  node  problem  is  seen  to  be  reduced  to  that  of  finding  that 

value  of  n  for  which  a  prescribed  function  F^  assumes  a  previously- 

computed  numerical,  value  .  This  problem  was  solved  by  inverse 

interpolation  between  known  values  of  F^  . 

13.  Once  nfl(t+6t)  ,  i.e.  n  at  the  node  g  at  time  t  +  fit  , 
p 

is  computed  the  cross-sectional  areas  of  the  stations  at  the  node,  i.e., 

A  „  ,  A.  ,  ,  and  A  may  be  computed,  and  then  finally  Q  n  , 
a, 2  d,1  c,l  a ,d 

a  ,  and  Q  ,  may  be  computed  from  the  approximate  continuity 
T>,1  c,l 

relations 
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(A28) 
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a, 2 


(t  +  St)  =  2Q 


a, 2 
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etc. 


(A29) 


lU.  It  is  important  to  note  that  the  mass  conservation  compatibil¬ 
ity  condition  expressed  by  Equation  A15  is  satisfied  by  the  Q’s  com¬ 
puted  from  Equation  A29.  An  interesting  property  of  the  computer  pro¬ 
gram  used  in  solving  the  channel  network  problem  is  that  any  of  the 
usual  errors  in  data  preparation  which  we  have  made  in  the  Savannah 
study  have  produced  numerical  outputs  in  which  Equation  A15  was  not 
satisfied  and  this  compatibility  condition  has  been  a  good  error  check. 


Network  Identification 

15.  Equations  A16-A27  were  given  for  a  specific  example,  i.e.,  the 
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typical  junction  indicated  in  Figure  Al.  In  application,  a  more 
general  procedure  has  been  adopted.  A  minimal  explanation  of  the  con¬ 
vention  used  in  network  identification  is  needed  to  facilitate  a 
description  of  this. 

1 6.  Every  station  was  given  a  unique  station  identification  such 
that  along  each  branch  the  station  ID's  formed  a  consecutive  sequence 
of  integers.  Branches  were  also  assigned  unique  integer  ID's  and  a 
table  of  the  first  and  last  stations  of  a  branch  versus  branch  ID  con¬ 
structed.  The  convention  that  the  positive  sense  of  a  branch  was  in  the 
direction  of  increasing  station  ID  was  adopted.  This  convention  fixes 
for  each  branch  the  direction  of  positive  flow  and  positive  slope. 
Finally,  a  table  of  branch  ID's  versus  node  ID  was  constructed  and 
was  used  with  the  convention  that  if  positive  displacement  along  a 
branch  implied  movement  toward  the  relevant  node  the  flux  into  the  node 
from  the  branch  was  positive.  The  sign  of  the  flux  for  the  various 
branches  was  also  tabulated  as  a  function  of  node  ID  . 

17.  Let  the  sign  of  the  node  flux  for  a  branch  a  at  node  $  be* 


aag 


+1  if  the  largest  station  number  for 
branch  a  is  at  node  B 

-1  if  the  lowest  station  number  for 
branch  a  is  at  node  B 


(A30) 


and  let  K  be  the  set  of  ID's  of  all  branches  ending  at  a  node.  Two 
other  sets  of  ID's  are  easily  determined.  These  are** 


M  =  {maB|m' 


the  maximum  station  ID  for  branch 


ending  at  node  B  if  =  1  ;  or  m^  =  the 

minimum  station  ID  if  =  -1^ 


(A31) 


*  The  o  table  must  be  constructed  prior  to  the  execution  of  the 
computer  program  and  is  input  (control)  data. 

**  The  long  |  is  known  as  the  universal  or  set  qualifier  and  may  be 
read  as  "such  that,"  the  {}  may  be  read  as  the  set  of  all  elements. 


(A32) 


N  =  <n  0  |n0=m0-a„> 
(  otB  1  aB  aB  aB| 


Using  Equations  A30 ,  A31,  and  A32,  the  generalizations  of  Equations  A22 
and  A23  are 


s»'  1 
neN 


S.  =  E  ^  o 
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where  neN  is  to  he  read  as  n  is  an  element  of  the  set  N 
tion  A25  is  generalized  as 


(A33) 
' (A3fc) 
Equa- 


F  (n)  =  I  «Am%o 
0  meM  m  2 


(A35) 


Equation  A26  is  unchanged,  i.e.. 


%  ■  S9  -  S* 


(A36) 


As  previously  explained,  the  equation  is  used  to  determine  n  and 
thence  the  A's  for  the  stations  at  the  nodes.  The  generalization  of 
Equation  A28  becomes 


6x 

-  ttt  (dA  +  6A  )  a  . 
2dt  m  n  aB 


(A3T) 


Equation  A29  is  unchanged,  i.e., 

^  (t  +  dt)  -  2^  -  Q^t)  (A38) 

18.  We  close  this  section  by  noting  that  the  essential  ideas  are 
presented  by  the  example  of  the  previous  section.  The  equations  pre¬ 
sented  here  are  given  as  a  matter  of  record. 

Boundary  Conditions 

19.  Various  approaches  to  the  treatment  of  boundary  conditions 


All 


were  tried,  but  ultimately  the  control  volume  approach  used  for  the 
junction  problem  was  adopted.  Assume  that  branch  a  has  an  extremity 
at  station  m  not  located  at  a  node.  We  define  o  by 


1  if  m  is  the  largest  ID  for  the  branch 


-1  if  m  is  the  lowest  ID  for  the  branch 


(A39) 


and  n  by 


n  =  m  -  a 


Using  the  approximation  of  Equation  A20  and  Equation  A17 


(6Am  +  6A  )  ~  =  cr  (Q  -  Q  ) 
m  n  2  a  n  m 


(AliO) 


(Ala) 


Station  m  is  at  the  boundary  and  either  n  or  Q  must  be  prescribed 
there  as  a  boundary  condition.  Also,  station  n  is  an  interior  station 
and  Aq  and  may  be  assumed  to  be  already  computed;  and  since  A 

is  a  prescribed  function  of  n  at  every  station,  either  or  Am 

is  also  known.  The  remaining  unknown  is  either 


6A  =  o  (Q  - 
m  o  n 


~  -  6A 
ox  n 

a 


(AU2) 


=  \  -  °o  (6Am  +  6An}  ~f 


(Ah3) 


where  Q^(t  +  6t)  may  be  computed  by  Equation  A29. 
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APPENDIX  B:  THE  SALINITY  MODEL 

Introduction 

1.  Numerous  problems  beset  the  predictive  determination  of  the 

salinity  distribution  in  a  real  channel  network.  In  particular,  the 
flow  is  turbulent  and  it  is  completely  infeasible  to  compute  with  a 
time  increment  <5t  and  space  increments,  <5x  ,  <5y  ,  6z  ,  which  would 

permit  the  solution  of  the  equations  that  describe  the  process  exactly. 

As  a  result,  simplifying  assumptions  and  correspondent  empirical  approxi¬ 
mations  must  be  introduced. 

2.  Techniques  for  accomplishing  the  necessary  simplification  and 
empiricism  have  been  the  subject  of  numerous  research  efforts.  A  brief 
review  of  the  mechanics  of  turbulent  diffusion  and  of  dispersive  pro¬ 
cesses  is  given  in  the  next  section.  It  is  hoped  that  this  discussion 
will  provide  some  insight  into  the  process  and  the  related  limitations 
of  the  resultant  mathematical  model. 

3.  The  one-dimensional  convective  dispersion  equation  is  usually 
used  to  describe  dispersive  processes  in  channels,  and  it  has  been  found 
that  for  a  proper  spatially  dependent  choice  of  the  dispersion  coeffi¬ 
cient  E  a  good  empirical  fit  to  observed  salinity  distributions  for 
various  estuarine  systems  may  be  obtained.  Unfortunately,  attempts  to 
correlate  the  dispersion  coefficient  with  channel  and  flow  characteris¬ 
tics  have  not  been  entirely  successful,  in  that  there  is  very  little 
coherence  between  estimated  dispersion  coefficients  and  proposed  empiri¬ 
cal  formulas.  Known  finite-difference  formulas  of  convective  processes 
suffer  from  the  phenomenon  of  numerical  dispersion.  Let  En  be  the 
numerical  dispersion  coefficient,  E  the  actual  dispersion  coefficient 

}  Q. 

which  exists  at  least  as  a  matter  of  formal  definition,  and  E  be  the 
dispersion  coefficient  arrived  at  by  empirical  data  analysis  or  mathe¬ 
matical  model  calibration.  Clearly, 


E  *  E  -  E 
e  a  n 


B1 


Also,  E^  is  dependent  on  the  numerical  model,  the  discretization 
parameters  6x  and  6t  used,  and  the  convective  rate  u  .  It  is  con¬ 
cluded  that  the  existence  of  and  variations  in  E  may  at  least  in  part 

n 

explain  variations  in  Eg  .  An  example  related  to  the  problem  of 
numerical  flood  routing  shoved  that  E  and  E  are  practically  equal. 

U.  The  great  empirical  variability  in  observed  disperion  coeffi¬ 
cients  raises  question”,  about  the  utility  of  numerical  predictions  of 
dispersion  processes.  Certainly,  results  must  be  interpreted  in  terms 
of  the  specific  problem  being  investigated  and  qualified  as  to  the  nature 
of  the  solution  provided.  In  this  context,  three  logical  questions  are 
significant : 

a.  How  sensitive  is  the  numerical  solution  to  variations  in 
the  dispersion  coefficient? 

b.  Do  numerical  solutions  indicate  the  same  relative  result  in¬ 
dependent  of  the  dispersion  coefficient? 

£.  Is  the  numerical  solution(s)  consistent  vith  the  behavior 
of  other  (hydraulic)  flov  characteristics  or  changes? 

The  extent  to  which  these  questions  may  be  affirmatively  answered  is 

indicative  of  the  reliability  that  may  be  associated  with  the  numerical 

solution. 


Turbulent  Diffusion  and  Dispersion  in  Channels 


5.  In  this  section,  an  idealized  channel  in  which  there  is  no 
lateral  velocity  is  assumed.  Let  x  be  the  horizontal  coordinate  and 
y  be  the  vertical  coordinate  and  u  and  v  be  the  velocities  in  these 
respective  directions.  Also,  let  s  be  the  concentration  of  a  solute 
which  is  conserved.  The  governing  equation  may  be  shown  to  be 


We  assume  that  we  may  decompose  s 


u  ,  and  v  into  component 


such  that  the  primed  part  is  that  which  is  due  to  turbulent  fluctua¬ 
tions.  In  the  manner  devised  by  Osborne  Reynolds  (see  Schlichting#) 
Equation  B1  may  then  be  transformed  into 


9s  9us  9vs  __  9u's 1  3v's 1 

9t  9x  9y  ~  “  9x  ~  9y 


(B3) 


For  practical  computation,  the  terms  on  the  right  of  Equation  B3  must 
be  empirically  approximated.  We  assume 


(Bit) 


and  the  diffusion  matrix  D  is  not  in  any  way  .dependent  on  s(x,y,t)  . 
Equation  BU  is  more  general  than  the  relations  which  are  usually  assumed, 
but  may  be  shown  to  be  less  general  than  the  actual  process,  i.e.. 
Equation  Bit  is  at  best  an  approximation. 

6.  It  will,  however,  suffice  in  our  investigation  of  the  main 
features  of  the  mechanics  of  the  diffusion  process  for  tidal  motion  in  a 
channel.  Let  us  now  nondimensionalize  Equations  B3  and  Bit  by  setting 


=  x*L  ,  y  =  y*h  ,  u  =  u^i/gh  ,  v  =  v*/gh  ,  t  =  t*L/»/gh  , 


Dij  =  d*,j  0  ® 


(B5) 


where  L  is  a  reference  length  and  h  is  a  reference  vertical  dimension 
and  the  scaling  parameter  o  is  for  now  left  undefined.  Combining 
Equations  B3  and  Bit  results  in 


9s  9us  3vs  9  /  9s 

9t  9x  9y  ~  9x  \D1,1  9x 


1,2 


Is  ’ \  +  JL  L  IS. 

3y/  9y  y  2,1  9x 


2,2 


(b6) 


where  the 


"  has  now  been  dropped. 


*  Dr.  Hermann  Schlichting.  1955.  Boundary  Layer  Theory,  Pergamon  Press. 
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g  Ei 
s  ii 


1.1  L  9x* 


+  °^D1,2  hip 


L  9y*  2,1  9x*  h  9y*  2,2  9y* 

where  e  =  h/L 

7.  For  tidal  flows  the  parameter  e  may  be  assumed  to  be  the 
depth  of  the  channel  divided  by  the  tidal  wavelength  and  thus  e  is  a 
very  small  number.  Also,  9s/9t*  ,  9s/9x*  ,  and  9s/9y*  should  be  the 

same  order  of  magnitude,  and  it  may  be  shown  that  ^  g  »  ^  i 

and  D  or  their  nondimensional  counterparts  are  of  the  same  order. 

2  *2 

An  examination  of  Equation  B8  then  leads  us  to  conclude  that 

a.  9v*s/9y  is  of  order  e  times  the  quantity  in  parentheses 
on  the  left  of  Equation  B8. 

b.  Turbulent  diffusion  is  of  considerable  importance  in  the 
dispersion  process  and  the  solution  of  Equation  B8  cannot 
be  in  reasonable  correspondence  to  the  physical  phenomena, 
unless  the  largest  turbulent  diffusion  term(s)  are  of  the 
same  order  of  magnitude  as  the  dominant  convective  terms(s) 

The  largest  diffusion  term  is 


£  9  __  9s 

h  9y*  U2,2  9y« 


and  observation  b,  therefore,  implies  that 


or 


o 


(BIO) 


The  resulting  first  order  model,  therefore,  is 


3s  3us  3_  n  _3s_ 
at  3x  *  3y  2,2  3y 


(Bll) 


8.  Pritchard*  has  made  an  extensive  analysis  of  field  data  ob¬ 
tained  from  the  James  River  estuary.  His  analysis  was  of  data  averaged 
over  a  tidal  cycle.  His  results  indicate  that  vertical  diffusion  was 
much  greater  than  horizontal  diffusion  and  that  horizontal  convection 
was  much  greater  than  vertical  convection  as  this  analysis  indicates 
they  should  be.  Li**  presents  an  analysis  of  a  number  of  not  basically 
dissimilar  problems  of  quasi-parallel  turbulent  flow  and  reaches  the 
same  conclusion. 

9.  To  obtain  a  better  feel  for  the  implications  of  Equation  Bll, 
let  us  first  assume  that  D2  2  equals  zero.  In  this  case. 


3s 

3t 


+ 


3us 

3x 


=  0 


and  the  behavior  of  s  is  as  shown  in  Figure  Bl.  In  particular,  for  a 
simple  shear  flow,  i.e.  u  dependent  only  on  y  ,  the  solution  may  be 
entirely  characterized  by  noting  that  along  the  curve 


x(y,  tQ  +  t)  =  u(y)  (tQ  +  x) 

s  as  a  function  of  y  is  exactly  the  same  as  the  distribution  of  s 
as  function  of  y  along  the  curve  x  =  0  at  time  t  . 


*  D.  W.  Pritchard.  1958.  "The  Equation  of  Mass  Continuity  and  Salt 
Continuity  in  Estuaries,"  Journal  of  Marine  Research,  Vol  17,  PP 
U12-U23. 

**  Li,  Wen-Hauing.  1972.  "Effects  of  Dispersion  on  DO  Sag  in  Uniform 
Plow,"  Proceedings,  American  Society  of  Civil  Engineers,  Vol  98, 

No.  SA1. 


B5 


*  =  o  at  time  t0- 


V"  u  =  o(y) 


Figure  Bl.  Graphical  solution  of  3s/3t  +  u  3s/3x  =  0  . 

For  a  simple  shear  flow 

10.  Secondly,  let  us  suppose  that  u  =  0  and  that  D_  0  is 

C. 

independent  of  y  .  Equation  Bll  reduces  to 

3s  _  afs 

3t  "  2,2  .  2 

3y 

and  since  there  can  be  no  flux  of  the  solute  across  the  bottom  of  the 
channel,  y  =  0  ,  or  through  the  free  surface  y  =  h 


—  I  -  IS.  |  =0 

*  U  ^  y=h 


A  solution  of  this  problem  is 


s  =  s.  + 
A 


^  e-DIjt  cos  x^y 


where 


a  1  n 

\) a  J  h 


Evidently  then  the  effect  of  vertical  diffusion  is  to  smooth  out 
whatever  variation  of  s  occurs  in  the  vertical  direction.  Horizontal 


convection  has  the  countering  effect  of  converting  longitudinal  gradients 
in  s  into  vertical  gradients. 

11.  Finally,  let  us  take  up  the  problem  of  obtaining  a  one- 
dimensional  model.  Integrating  Equation  Bll  with  respect  to  y  results 
in 


3s  3_  1 
3t  3x  h 


usdy 


D2,2  3s 

h  3y 


y=h 


y=0 


(B12) 


where  the  superscript  ~  denotes  an  average  over  the  depth.  The  no¬ 
flux  condition  at  the  free  surface  and  bottom  of  the  channel  reduces 
the  right-hand  side  of  Equation  B12  to  zero,  or 

,» 

I t'ktj  *•*-»  <B13> 

o 

As  a  matter  of  formal  definition.  Equation  B13  may  be  written  as  (as  it 
has  been  in  practical  applications) 


3s 

3t 


+  u 


3s  _  3  3s 

3x  3x  3x 


(BlU) 


Equation  BlU  is  disturbing,  because  there  is  no  analytical  technique 
for  obtaining  it  from  Equation  B13  and,  therefore,  no  theoretical  reason 
to  anticipate  that  E  is  a  function  of  some  set  of  variables.  Equa¬ 
tion  B12  is  a  diffusion  equation,  whereas  Equation  BlU  is  a  dispersion 
equation.  It  should  be  noted,  however,  that  the  combination  of  a  verti¬ 
cal  velocity  gradient  (i.e.,  a  shear  flow)  and  vertical  diffusion  cause 
an  apparent  horizontal  dispersion.  Since  no  theoretical  definition  of 
E  is  possible  it  may  be  viewed  as  an  "empirical,  knob"  which  when  prop¬ 
erly  twisted  will  bring  the  solution  of  Equation  BlU  into  a  best- 
fit  correspondence  with  a  given  set  or  sets  of  field  data.  Nothing  of 
course  is  wrong  with  this  if  a  proper  prescription  for  twisting  the  knob 
can  be  found  and  if  the  approximation  is  not  too  sensitive  to  minor 
mlsadjustment . 


BT 


A  Numerical  Approximation 

12.  The  one-dimensional  model  of  dispersion  for  a  conservative 
constituent  s  ,  which  will  now  be  identified  as  the  saline  concentra¬ 
tion,  is 


9s  3s  9  9s 

9t  U  3x  9x  9x 


(B15) 


13.  Assuming  that  E  is  zero,  an  approximate  solution  of 
Equation  B15  is 


s(xi,  tJ+1)  =  s(xi  -  uisj+1  tj)  (Bl6) 

In  the  numerical  computation  u  and  s  are  known  only  at  the  stations 
along  a  branch.  Let 

sij  =  s(v  V  “d  "i,j  ■  u(v  V  (B17) 

6t  was  chosen  so  that 

u6t  <  6x  (Bl8) 

for  all  branches;  and  s.  .  .  was  found  by  a  linear  interpolation 

1  » 

between  s.  ,  and  s.  .  if  u.  .  .  was  positive  or  s.  .  and 

lU.  The  algorithm  Just  indicated  is  subject  to  numerical,  dispersion, 
and  it  was  not  found  to  be  essential  to  add  "real"  dispersion  in  order 
to  obtain  good  agreement  between  the  model  and  prototype  for  the 
Savannah  system.  In  the  Savannah  study,  the  model  was  used  to  compute 
steady-state  salinity  conditions  (i.e.,  repetitive  tide  and  salinity 
at  boundary)  rather  than  transient  conditions. 

Junctions  and  Boundaries 
15.  The  one-dimensional  dispersion  equation  was 


B8 


I 


3s  .  3s  _  3  _  3s 
3t  U  3x  "  3x  E  3x 


(B19) 


while  the  mass  conservation  equation  of  an  incompressible  fluid  is 


o 

3t  3x 


(B20) 


multiplying  Equation  B19  by  A  and  Equation  B20  by  s  and  adding 
results  in 


(A  It  +  S  It )  +  (Q  If  +  S  If)  =  A  fe 


E 


3s 

3x 


or 


Ms  +  3Ss  =  Ai-Ei£ 

3t  3x  3x  3x 


If  dispersion  is  ignored  (i.e.,  E  =  0),  Equation  B22  reduces  to 


Ms  3$s 
3t  3x 


(B21) 


(B22) 


(B23) 


which  is  essentially  the  same  as  Equation  B20. 

1 6.  Appropriate  compatibility  conditions  at  a  node  are 


Qs  (from  all  branches)  =  0 


and 


(B2U) 


s(branch  1  node  station)  =  s(branch  2 —  n  node  stations)  (B25) 

Recalling  that  the  compatible  conditions  used  in  determining  the  flow  at 
a  Junction  (pages  A6-A11)  are  essentially  the  same,  one  sees  that  the 
same  algorith  may  be  applied. 
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One  obtains 


Or,  since  6Fn  is  in  this  case  known 


where  like  and  S„  ,  S'  and  S'  have  been  previously  computed 

fed  H  *r* 

17.  The  same  correspondence  exists  at  boundaries,  so  they  will 
not  be  discussed  separately. 
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